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Abstract 

We propose a straightforward and model independent methodology for characterizing the sensi- 
tivity of CMB and other experiments to wiggles, irregularities, and features in the primordial power 
spectrum. Assuming that the primordial cosmological perturbations are adiabatic, we present a 
function space generalization of the usual Fisher matrix formalism, applied to a CMB experiment 
resembling Planck with and without ancillary data. This work is closely related to other work 
on recovering the inflationary potential and exploring specific models of non-minimal, or perhaps 
baroque, primordial power spectra. The approach adopted here, however, most directly expresses 
what the data is really telling us. We explore in detail the structure of the available information 
and quantify exactly what features can be reconstructed and at what statistical significance. 



I. INTRODUCTION 



In the last decade, CMB experiments have created new opportunities for studying the 
early universe. With the data obtained by WMAP [lj and the promise of even more accurate 
CMB measurements by Planck j2] in addition to the several current and future ground and 
balloon based experiments, the data will become even more constraining and it will be 
interesting to see whether the simple minimal model that seems to be able to explain all the 
present data survives or whether a new layer of complexity is discovered. 

According to our current best understanding of primordial cosmology, the universe un- 
derwent a period of extremely rapid expansion known as inflation [3H6]. This model was 
initially proposed to solve several problems with the standard big bang theory. In addition, 
inflation was found to make quantitative predictions concerning the primordial cosmological 
perturbations [7HIU]. In the simplest of these models, inflation predicts that the primordial 
fluctuations in primordial power spectrum (PPS) should be Gaussian with a nearly scale 
invariant spectrum. The most common method of interpreting the CMB data assumes a 
modest number of cosmological parameters (for example, in the WMAP 7-year analysis the 
parameters adopted are Hq, Ub, uj c , Q\) and a parameterized form for the primordial power 
spectrum, for example P(k) = A(k/k ) ris ^ 1 or perhaps P(k) = A{kj 'k p n3 ~ 1 ~P ln ( k ' k °>> where 
the parameter (3 denotes the running of the spectral index. In this paper we adopt as a fidu- 
cial reference model the WMAP 7-year best fit [T], for which the parameters are h = 0.703, 
uj h = h 2 Q b = 0.02227, u c = h 2 Q c = 0.1116, r = 0.085, A = 2.42 x 10~ 9 , and n s = 0.966. 
Here the somewhat arbitrary pivot scale has been set to ko = 2.0 x 10 -3 Mpc~ l . Further- 
more we assume that there are no tensor perturbations nor running in the scalar spectral 
index (i.e. r = and /3 = 0.0) and also that the universe is exactly flat. It is encouraging 
to learn that a simple model of this sort can provide an adequate and plausible explanation 
of the presently available data, and following this general approach, one can construct more 
complicated models and compare these based on relative likelihood and ideas about model 
comparison. 

However the drawback of an analysis comparing parameterized models to the data and 
among themselves is that which aspects of the model actually enter into the likelihood 
remains hidden inside a black box. If for example we show plots of constraints on n s or 
the running /3, these plots do not show at what k measurements have been made. The 
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purpose of this paper is to enter into this black box and show numerically what exactly is 
being measured, in particular with respect to the primordial power spectrum. For simplicity 
we assume that the perturbations are adiabatic and Gaussian, although generalizing the 
treatment here to isocurvature perturbations is in principle straightforward. 

Rather than developing a suite of fancy parametric models with lots of bells and whistles, 
we approach the problem more in the spirit of the method dubbed 'exploratory data analysis' 
by John Tukey (see for example [11]) where the aim is to devise methods to make manifest 
patterns in the data and what is being constrained and where, the purpose being to uncover 
the unexpected. This approach stands in contrast to the complementary approach where 
precise models with well-defined prior information are set forth as hypotheses and optimal 
statistical tests to discriminate between them are sought. In fact the use of priors in the paper 
is more formal and metaphorical, because for our purposes meaningful prior information is 
lacking. The priors serve more as regulators to suppress noise where there is seemingly no 
relevant information. 

In the past a number of approaches have been proposed to search for features in the 

primordial power spectrum inferred from the CMB. One general class of approaches consists 

of parametric models, where P(k) is modeled using a family of functions having a fixed and 

well-defined number of adjustable parameters. Bridle et al. [T2], for example, analyzed the 

WMAP first-year results using a model with N knots equally spaced in ln(fc) between a k m i n 

and k max where \n(P(k)) is interpolated linearly between N control points or knots. The 

number of knots admitted is a measure of the model complexity, and as is well-known, when 

comparing models it is necessary to compensate for the better fit due to fitting away noise 

when more degrees of freedom are present. This can be done using standard techniques (e.g., 

a simple x 2 analysis, the Aikake information criterion (AIC) [13j, the Bayesian information 

criterion (BIC) |14^] , etc). One disadvantage of this approach is that the placement of 

the bin boundaries, which is arbitrary, influences the conclusions. The sensitivity to a sharp 

feature having a width comparable to the bin width depends on where the feature is situated 

relative to the bin boundaries. Interesting variants have been proposed. Hamann et al. have 

proposed allowing the knot positions not to be fixed but rather parameters of the model 

1 Care must be taken when applying the BIC to actual data because of the asymptotic nature of its 
derivation. It is often asked what value to use for ndata in the log(rid ata ) term appearing in the BIC. Is 
it the number of Ci or the total number of measurements? A careful analysis starting from a Bayesian 
prior reveals that while the BIC obtains the right scaling as ndata ~~ ^ oOj there is another prefactor that 
does not vary with ndata that is ignored. Without knowing this prefactor, which would required a specific 
Bayesian prior, it is not possible to address this question. In reality, the factor should be ln(rod ata /n) 
where n depends on the overall relative likelihood of the models in the prior. 
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[15]. Although this approach introduces nonlinearity into the model, it allows localized sharp 
features to be modeled by placing control points only where they are needed. In this way 
models giving a good fit to a single sharp feature are significantly less penalized than in the 
approach of Bridle et al. and others [T2l IT6] . 

Another class of approaches may be described as penalized likelihood methods. In this 
class of methods the function space for P(k) includes many more parameters than desired. 
There is either a very large or infinite number of degrees of freedom. Model complexity is 
limited by penalizing very rough or rapidly varying functions. The log likelihood is modified 
to include a roughness penalty for example as in the following form: 



Here we may appeal to Bayesian doctrine by describing the second term as a prior probability 
over the function space of P(k). An example of this approach applied to the WMAP data 
can be found in [17]. However in practice this term rarely corresponds to any actual prior 
belief. Rather it is chosen a posteriori to obtain a meaningful result and depends on the 
quality of the data. The form of the roughness penalty in eqn. ([T| is somewhat arbitrary 
and other reasonable forms have been proposed. We could also have used the square of 
the first rather than the second derivative, although this choice would have the drawback 
of slightly favoring an overall flatter spectrum. Also the weighting as a function of ln(fc) 
and what is assumed for k smaller and larger than where there is data involves arbitrary 
or subjective choices. The challenge of choosing the regularization parameter A involves 
finding the right compromise between underfitting and overfitting the data. There exists an 
abundant literature on how to make an optimal choice for A rapidly surveyed below. However 
in our view it does not suffice to find the optimal A, because this information does not tell us 
which features of the reconstruction are to be believed and which are to be regarded as most 
likely explained as statistical fluctuations. The approach adopted in the paper is to study 
and characterize precisely the properties of the reconstruction for a variety of reasonable 
values for the smoothing parameter. 

Tocchini-Valentini et al. [TBI [19] adopted a penalized likelihood similar to that used 
in this paper, which was applied to the WMAP data finding some evidence for statistically 
significant features. In this paper we account for uncertainties in the cosmological parameters 
rather than fixing them and adopt a frequentist approach to the error analysis because we 




(1) 
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do not believe that the prior implicit in the smoothing penalty deserves to be interpreted as 
an actual prior belief or probability measure on the function space of possible power spectra. 

Many reconstruction procedures can in some cases yield reconstructed primordial power 
spectra that take negative values in some places, which is clearly unphysical. This can occur 
when a quadratic approximation to the likelihood is made or when there is a superimposed 
noise component. Hamann et al. [15] (see also [20]) carried out a modified Lucy- Richardson 
deconvolution of the primordial power spectrum from CMB data. Lucy [21] and Richardson 
[22] introduced an algorithm widely used for image deconvolution having the property that 
the reconstruction is always positive. The original Lucy-Richardson algorithm was derived 
from a likelihood assuming a Poisson noise for the number counts. If there are some counts 
in a cell, an infinitely improbable likelihood results as the predicted average count tends to 
zero, quite unlike in a Gaussian, where a negative average count would not be excluded by 
the form of the likelihood. Unfortunately the CMB likelihood does not fit into this mold. 
Modifications to L-R algorithm have been proposed to take into account the properties of the 
CMB likelihood. However the modified L-R algorithm cannot be derived rigorously starting 
from a CMB likelihood. Hamann et al. overcome this difficulty by carrying out a frequentist 
analysis of the reconstructed power spectra, treating the reconstruction as a black box and 
characterizing its properties through Monte Carlo experiments. 

Hu [23] and Leach [24] explore an approach based on principal component analysis. A 
principal component analysis ranks directions of deviation with respect to a fiducial power 
spectrum using the eigenvalues and eigenvectors of the Fisher information matrix. The 
decomposition into eigenmodes, however, relies on a notion of length (or metric) for this space 
of linearized deviations, and thus entails some arbitrariness]^] Peiris et al. [T7J[25ll26] explore 
a penalized likelihood using 'cross-validation' to select an optimal degree of smoothing. A 
variety of other approaches have been proposed including wavelets (Mukherjee et al. [271429] ) 
and 'cosmic inversion' (Kogo et al. [3D], Nagata et al. [3H |32] and Ichiki et al. [33]). 

Vazquez et al. [31] present an interesting analysis based on Bayesian model comparison. 
A prior distribution is constructed for a linear interpolation of ln(P) as a function of ln(/c) 
using a variable number of nodes whose positions are free, yielding a posterior distribution 
among the models. Using data from WMAP combined with ACT and SDSS LRG data, the 
2 A quadratic form is a linear mapping from a linear space into its dual and not into itself. Therefore a 
metric tensor is required to define a natural equivalence between the linear space and its dual. 
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authors conclude that an ns = 1 Harrison- Zeldovich- Peebles spectrum is 'strongly' excluded 
within a Bayesian framework, using the terminology of Jeffreys. 

This paper systematically studies and characterizes the trade-offs between bias and vari- 
ance (i.e., noise) in the choice of smoothing scheme. Rather than espousing a particular 
philosophy leading to a unique solution, we show in a quantitative way the results for a 
range of choices. Which is best depends on what type of features one is trying to recon- 
struct. 

The organization of the paper is as follows. Section 2 presents a detailed description of the 
constraints on the primordial power spectrum assuming a Planck-like experiment based on 
the Fisher information kernel. The information loss due to uncertainty in the determination 
of the other cosmological parameters is described. Section 3 describes the power spectrum 
construction first for the simplified case where the other cosmological parameters are fixed 
and then including their uncertainties. Section 4 characterizes the noise of the reconstruction 
and with what statistical significance various model features can be reconstructed. Finally 
Section 5 presents some concluding remarks. 
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II. FISHER INFORMATION KERNEL 
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FIG. 1. Fisher information density relative to a reference primordial power spectrum. 

The information density defined in eqn. ([6]) is shown. 



The full-sky temperature power spectrum Ci depends linearly on the primordial power 
spectrum and is defined as the following linear integral transform 

/oo 
d\nkP(k)[A ( p(k)} 2 (2) 
■oo 

where A) (k) is an adiabatic temperature transfer function as computed using standard 
Boltzmann solver such as CAMB [35J. The polarization and cross spectra, Cf E and Cj E , 
are similarly defined. We approximate the likelihood by assuming full sky coverage and 
uncorrelated white noise 
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(obs) 



(obs),AB 



}a,b=t,e is the measured power spectrum matrix, the the- 



oretical (predicted) power spectrum matrix, and is the detector noise. In this section 



j-1 I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I 

10" 4 10" 3 10" 2 10" 1 

k {M P c~ l ) 




k c (MpcT 1 ) 

FIG. 2. Properties of the Fisher information density kernel. The top plot shows the 
width as defined in eqn. Q. The middle plot shows F(k) as defined in eqn. J9]), which provides an 
estimate of how accurately fluctuations at the scale k are being constrained. The bottom plot shows 
the structure of the correlation matrix (i.e., the kernel C(k, k') = I(k, k')/[I(k, k)I(k'k')] 1 / 2 ). Here 

horizontal coordinate is fc^ or = s/kk' and the vertical coordinate is \n(k'/k). On the far left and 
the far right the elements far away from the diagonal have a correlation of almost one, indicating 
that only a single quantity is being measured at very small and very large k. 



and elsewhere in the paper unless otherwise specified, we assume a fictitious experiment 
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somewhat similar to Planck based on the specifications published in the bluebook [2] rather 
than on an updated understanding of the mission and the instrument in-flight performance. 

In this paper we are concerned with the effect that small scale variations in the primordial 
power spectrum have on the likelihood. We define the Fisher information kernel I as 

-21nL = J dink J d\nk' [P{k) - P ML {k)}I{k,k')[P{k') - P ML {k')} (4) 

where an irrelevant constant term has been omitted. Using this definition we can express 
the kernel as a functional derivative of the log likelihood 



/(In k, In k') 



5\-\nL) 



(5) 

P=Pml 



SP(\nk)SP(\n k') 

where Pml is the primordial power spectrum at maximum likelihood. 

Since we already have a rough idea of what the primordial power spectrum looks like 
on the scales probed by the CMB, it makes sense to describe the power spectra under 
consideration in terms of fractional variations with respect to a fiducial power law spectrum 
(for which the parameters are taken from the WMAP 7-year best fit model). We set Po(k) = 
A(k/k ) ns ^ 1 where A = 2.42 x 10~ 9 , k = 2.0 x 10~ 3 Mpc' 1 , and n s = 0.966, and define 
P{k) = Po(k)[l + f{k)]. Let us for the moment fix the other cosmological parameters to 
their WMAP 7-year best fit values. We may express the expectation value of the relative 
log likelihood as the following quadratic form about the best-fit value as follows: 

A[-lnL] = l -J d(lnh) J d{\nk 2 )[f(h) - f M L(h)]I(h,k 2 )[f(k 2 ) - f ML (k 2 )] (6) 

where I(k,k') = P (k')I(k, k')P (k). In our definitions we have tried to remove dimensions 
from quantities wherever possible. Fig. [l] shows the Fisher information density I(ki,k 2 ). 
The integral over both arguments is interpreted as the \ 2 f° r detecting a CMB signal over 
a null hypothesis of no CMB signal at all, and as expected this quantity is approximately 
10 6 or t 2 max where £ max is approximately the multipole number where the signal-to-noise 
ratio {S/N) has dropped to one. Were it not for correlations in the Fisher information 
density, we would be able to measure about a million independent quantities characterizing 
the shape of the power spectrum at a marginal statistical significance — that is, at about 
la. This would indeed be the case if the Fisher information density were more strongly 
peaked along the diagonal. The width of the peaked region along the diagonal provides an 
approximate idea of how coarsely one has to bin for correlations between adjacent bins to be 



negligible. If the peaked region along the diagonal is very wide, we measure a small number 
of quantities (much smaller than Xtotai) & lt>eit at very high precision. The most extreme 
case would be a Fisher matrix where all the entries are equal, or more generally of the 
form J(ki, K2) = g{K>\)g(K>2) where we define k = ln(k). In this case the Fisher information 
operator would be of rank one and precisely one quantity can be being measured. Fig. [2] 
shows the diagonal of the Fisher information kernel and the correlation matrix obtained from 
the Fisher information kernel — that is C(k, k') = I(k, k' )/[I(k,k)I(k , ,k')} 1 / 2 . The fact that 
at very low k and very high k the correlation becomes very nearly one even quite far from 
the diagonal indicates that in each of these two region there is only one relevant quantity 
that is being measured constituting a sort of weighted average of the power spectrum at 
each end. It is convenient to define the local width of the diagonal in the following way 
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FIG. 3. Fisher density kernel transverse profile. We plot the Fisher density kernel in 
the direction orthogonal to the diagonal (plots labeled according to wavenumber k c of diagonal 

crossing), (i.e., I(k c + AA;, k c — Ak) 
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to (In k) 



Jd(\nk') I(lnJfe,lnJfe') 



(7) 



/(In k, In k) 

This quantity is plotted in Fig. [2] in the middle panel. For A(ln/c) ;> w, the error in the 
binned power spectrum may be approximated as 

~5P 1 1 \ (\ \ /i-M-Va 

a — « / -lnfcifc 2 , -lnfcifc 2 Uu -lnfci/ 

_ -r jfce[fei,jfc 2 ]J |_ \2 2 / \2 

Fig. [2] shows (on a logarithmic scale) 

-1/2 



(8) 



F(ln fc) = /(In fc, In k)w(\n k) 



(9) 



When / ^ 1 we have essentially no useful information regarding the power spectrum given 
current expectations. 
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FIG. 4. For the partition into bins indicated at the bottom, the images show the correlations 
between the bins using a greyscale. 100 bins were chosen so that the information in each bin as 
defined in the text is approximately equal. Each bin contains a x 2 corresponding to a detection at 
about IOOct of power within that bin. 



We may choose bins of varying width so that each bin contains roughly the same amount of 
information. The bottom of Fig. [4] shows how the domain in wavenumber can be subdivided 
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into approximately 100 bins each containing about a detection at 100a (relative to the null 
hypothesis of a vanishing power spectrum in that bin). The locations of the bin boundaries 
are shown in the lower panel of the plot. The upper panel of Fig. [4] shows using a greyscale 
density plot the correlations between bins. 

On smaller scales the story is more complicated. The precise profile of the correlations 
near the diagonal is relevant for determining what information can be extracted concerning 
small-scale structure in the power spectrum. In Fig. [3] the transverse profile of the Fisher 
information density kernel is shown for some representative wavenumbers. The shape of the 
transverse profile is key to determining how much information is available for reconstructing 
variations in the primordial power spectrum that change rapidly with k. A smooth broad 
profile acts to erase information on small scales in k whereas a sharp narrow profile facilitates 
recovering information concerning such rapid variations. 

Relative to a given experiment or collection of datasets, over what range in k are the 
parameters of the power spectrum, here 11$ and /3, actually being measured? This is an 
important question because the parametric model serves more as a local fitting function 
for the power spectrum and does not derive from a fundamental theory allowing us to 
shamelessly extrapolate the power law as far as we want without trepidation. If we really 
knew that our power law ansatz were correct, this would not be a relevant question. We 
would like to define k 1ow and K hl9h to characterize how far to the left and to the right the bulk 
of the statistical information extends. Given the parameterization for the power spectrum 

P(k) = A(k/k ) {ns+l3Hk/ko) \ (10) 

it follows that 

\nP(k) = \nA + n s \n(k/k ) + f3\n 2 (k/k ) (11) 

and to lowest order 

V = V + (Sn s )(K - «o) + (*/?)(« - K Q f (12) 

where k = Ink, = lnfco, and ko is the so-called pivot or reference scale. Under a change 
of pivot scale Kq — > Kq, the parameters transform as 

A(k ) -> A(k ) = exp [lnA(/c ) + n s (K )(k - k ) + P(kq)(k - k ) 2 ] , 
nsM -> ^s(£o) = n,gf(«o) + /9(«o)(«o - K o), 

12 



£(«„) P(ko) = /9(«o) (13) 

where Ko — lnfco and Ko = ln/co- It is sensible to choose ko so that the uncertainties in A 
and ns are uncorrelated. Otherwise, because of a lever arm effect, the errors become large 
and highly correlated. Alternatively, the lack of correlation between ns and /3 could have 
been chosen to set Kq, but we have only one free parameter to adjust, and both conditions 
place ko within what may be described as the 'center of mass' of the statistical information 
concerning the form of the primordial power spectrum for a specific experiment or a collection 
of datasets combining multiple experiments. 

Having chosen a sensible pivot scale, we first consider a simplified situation with a diagonal 
Fisher information density for the power spectrum so that 

X 2 = J dK Idiag(n)[Pmodel(n) ~ Pdata{^)f ■ (14) 

In this situation, the condition fixing the pivot scale is 

j tk Idiag{K){K - Kq) = (15) 

and the unmarginalized information (inverse of the variance of the measurement with the 
other parameters fixed) ns and (5 are given by the two integrals 

In S = J dK I diag {K){K - Kq) 2 (16) 

and 

1(3 = J dK I diag {K){K ~ Kq) 4 . (17) 

For the Fisher density of an actual experiment, there are always correlations even if the 



information is strongly clustered near the diagonal and the integrals in eqns. (16) and (17) 



become double integrals. We modify eqns. (16) and (17) by limiting the integration over 



just one of the two variables of integration. For the information concerning each variable, 
the cumulative information up to a certain k is defined as follows, with the expression for 
ns given by the integral 

/k /*+oo 
cLk I d,K I(k — K /2, K + k /2)(k — Ko) 2 , (18) 
-oo J — oo 

and we define k 1 °™ and k^I h as the locations of the the first and third quartiles, respectively, 

so that In s {^ns) = (V 4 )-^ and J» 5 («n/ fc ) = ( 3 / 4 ) J « S ; and 4™ and ^ h can be defined 
analogously. 
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The shape of the cumulative information for ns and are shown in Fig. [5] as well as 
the locations of the quartiles. (Note that when the matrix is not concentrated along the 
diagonal and especially in the presence of negative correlations assigning information to a 
specific k no longer makes sense and the above definitions become problematic.) 





FIG. 5. Cumulative information for ns and (3. The left panel shows the cumulative information 
for ns as defined in eqn. ( 18 ) as well as the first, second, and third quartiles. The right panel shows 
the analogous information tor the running in the spectral index (3. 
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III. POWER SPECTRUM RECONSTRUCTION 



A. Basic approach 

The previous section described the Fisher information density for the fractional variation 

/(*) = 6P(k)/P (k) (19) 

of the primordial power spectrum relative to a reference model taken to be the 7-year 
WMAP best-fit model p] already described in the previous section. The Fisher information 
density I(k, k') was considered both with the vector (3 containing the other cosmological 
parameters fixed — that is, assumed to be determined extremely precisely from other non- 
CMB datasets — and with (3 determined internally using the CMB data alone and hence 
uncertain. In the latter case we marginalize over the other cosmological parameters treating 
them as nuisance parameters with a flat (i.e., non-informative) prior probability distribution. 
In this section we consider how best to reconstruct the power spectrum — or more precisely 
f(k) — over the range of wavenumbers where the CMB provides sufficient information. 

In principle, if the parameters (or perhaps bins) used to describe f(k) are less numerous 
than the number of independent data points, we could reconstruct f(k) without a regulator, 
or prior distribution concerning the smoothness of f(k) in the form of a roughness penalty. 
However since we want to avoid a parametric model as well as artifacts from bin boundaries, 
we instead introduce a roughness penalty. The Fisher information matrix for the primordial 
power spectrum takes the general form 

/,/ = (X T ) fc I cc X cf + XRff. (20) 

In the quadratic approximation, minus twice the log likelihood is given by x 2 — f T -^//f • Here 
I cc is Fisher matrix for the angular CMB power spectrum. In order to provide as compact a 
notation as possible, we write our expressions as if we were dealing with finite-dimensional 
matrices and vectors, even though the index / represents continuous variable and c is a 
discrete index. Here XRff represents a regulator whose eigenvalues should very nearly 
vanish for slow variations in the power spectrum but be large (relative to the eigenvalues 
of the Fisher matrix) for modes where f(k) wiggles rapidly. A commonly used but by no 
means unique choice for this regulator is the quadratic form 

XR(f) = Xi T R f f{ = xfdK (^^) 2 , (21) 
15 



which is a measure of how much the function f(k) deviates from a straight line. Under this 
regulator, to linear order a change in the amplitude A s or spectral tilt rig suffers no penalty 
at all, but rapid variations in the running of the spectral index, or a large curvature, are 
penalized by an amount proportional to the metaparameter A. The object is to penalize 
functions that are too rough through a judicious choice for A. 

Above we have assumed that the other cosmological parameters, or elements of the vector 
/3, have been fixed. If instead these are determined internally using the CMB data, the ex- 
pression for the Fisher information for the primordial power spectrum after marginalization 
is modified to 

j-margjlat _ (X T ) fc [l cc - (I cc X c/3 ) (Xj c / cc X c/3 ) _1 (X^J CC )] X cf , (22) 

and in the presence of additional information on (3 (for example, by incorporating other 
external non-CMB data sets) expressed by the Fisher matrix Jg^, this expression modifies 
to 

j(rnarg, prior) _ (X T ) f c [/ cc - (I cc X c/3 ) (Xj c I cc X cf3 + Ipp) _1 (X^J CC )] X c f. (23) 

The equations here assume a Gaussian likelihood as well as Gaussian prior information. The 
most significant non-Gaussianity arises from the angular power spectrum likelihood at very 
low i and from nonlinearity in the change of the angular power spectrum as a function of 
variations in the cosmological parameters (i.e., the components of the vector (3). Supposing 
Gaussianity allows simple analytic expressions to be obtained. In practice the distributions 
are slightly non-Gaussian and must be explored using MCMC or some other numerical 
sampling method. The Gaussian analysis, however, is invaluable for gaining intuition about 
how to select a good regulator, but the solutions obtained should subsequently be validated 
by simulations taking into account the non-Gaussianity of the regulated likelihood. 

Fig. [6] shows the Fisher density when h, Ub, and u c are uncertain and constrained only by 
the CMB data, and Fig. [7] shows some transverse profiles. These should be compared with 
the analogous plots in Figs. [T] and [3] of the unmarginalized Fisher density, which are almost 
identical at small k, but have significantly different behavior at high k. Marginalization 
over the cosmological parameters increases the inter-A; correlations, as seen in Fig. [8] In 
marginalizing over the cosmological parameters, we have not included the reionization optical 
depth r, which is doubtless a relevant parameter. This is because except at very low £, 
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FIG. 6. Fisher density contour plot (after marginalization over the cosmological parameters). 

the reionization optical depth r and the overall amplitude of the power spectrum A s are 
completely degenerate. It is the combination A s e~ 2r that normalizes the overall amplitude 
but not the shape of the CMB sky power spectrum. Fig. [9] shows the deviations to this 
statement at low t. The deviation is greatest at low £, especially for the E-mode polarization. 
The data from the very first multipoles can be used to determine r, and then to determine 
the overall normalization. 

We now return to the problem of choosing the metaparameter A of the regulator. On 
the one hand, if A is too small, the reconstruction is dominated by noise, seeded by the 
cosmic variance of a specific sky realization and by random instrument noise on smaller 
scales. In this situation, the bias of the reconstruction will be small but its variance will be 
large. On the other hand, if A is too large, real features in the primordial power spectrum, 
especially those that are sharp or highly localized, will be smoothed over and distorted. 
This is bias accompanied by very little variance. Finding the right A amounts to making the 
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FIG. 7. Cross sections of the Fisher information density along transverse direction 
(with cosmological parameters uncertain and determined from CMB data alone) cen- 
tered on the diagonal. The k c indicate the diagonal crossings. Compare these plots to Fig. [3j At 
large scales (log 10 k c = —4 and log 10 k c = —3) the kernel is virtually identical in both figures, while 
at small scales the behavior changes depending on whether the cosmological parameters are fixed 
or internally determined with uncertainty. 

right compromise between bias and variance, both of which should be minimized, resulting 
in conflicting requirements that cannot be satisfied simultaneously. 

The effect of the regulator may be described using the language of Bayesian statistics, 
as a prior distribution over a function space. While this description provides useful insight, 
the resulting smoothness prior does not correspond to any actual prior belief concerning the 
theoretical power spectra. Rather the prior is chosen after the fact and is informed primarily 
by the characteristics of the experiment and the information provided by the data. 

Automatic procedures can be derived for choosing the right parameters for the regulator. 
In the course of this study we tried using the Aikake Information Criterion (AIC) p2] 
to choose A as well as criteria penalizing model complexity more strongly inspired by the 
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FIG. 8. Correlation after marginalization over the uncertain cosmological parameters 
determined from CMB alone. Note that there is now much more correlation between very low 
and very high k. 
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FIG. 9. Ratio of the temperature and E-polarization C-g's at two different values of the optical 
depth: t\ = 0.1, T2 = 0.08. For t <; 30, the shapes differ by less than 1%, and for t J> 100, less 
than 0.1%. 



BIC (Bayesian Information Criterion) [Hj. However, while these criteria may offer some 
theoretical rational for selecting A, they do not provide any insight into what kind of features 
in the power spectrum can and cannot be recovered. 
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The problem of choosing a regulator to recover the primordial power spectrum necessarily 
involves some degree of arbitrariness. The approach adopted here is pragmatic. We construct 
a family of model features to be reconstructed and then test various regulators with a range 
of smoothing parameter choices. The final product puts forth a reconstruction procedure 
whose properties have been well characterized by showing what features can and cannot be 
recovered given the experiment and at what statistical significance. 



The functional form of the regulator adopted in eqn. (21) is somewhat arbitrary. This 
form could be justified if there were a translation invariance in the variable k = Ink. The 
discussion of the previous section showed that the CMB provides stringent constraints on 
f(k) over a window spanning approximately two decades in k, but outside this window 
f(k) is hardly constrained. In fact, for k well below this window — that is, for modes whose 
wavelengths are superhorizon today — we measure only the local curvature in the geometry 
and have no way of deducing the precise wavelengths of the long-wavelength modes giving rise 
to this curvature, which may be approximated within the present horizon using a quadratic 
approximation. This situation implies that only an integral over the very \ow-£ power 
spectrum is constrained by the CMB data. This fact is evident from the correlation structure 
of the Fisher matrix at low-fc (Figs. [2j [8]), which shows that the dimensionless correlation is 
almost one in the corner. Because of the sharpness of the centrifugal barrier in the differential 
equation for the spherical Bessel functions, the very \ow-k power spectrum affects only the 
quadrupole and has almost no impact on the higher-^ multipoles. The impact of the high-£ 
part of the power spectrum lying above this window is more complicated because at large 
arguments the spherical Bessel functions oscillate rapidly and fall off slowly according to a 
power law that does not depend on £. 

This situation suggests that more regulation is needed outside this window than inside, 
suggesting the following generalization of eqn. pi) 

}_V(^ 2 
dn 2 

yet this proposal introduces additional arbitrariness in the choice of p. In order to deal with 
the region of the power spectrum that is hardly constrained by the CMB data, we found 
it necessary to modify the regulator by adding a term pushing the reconstructed spectrum 

toward the fiducial model. Thus the regulator becomes 

r /8 2 f(K)\ 2 f Kmin r +co 
f T Rjf(\, a)f = A / dK ( ^ 2 ' J +a dnf^ + a ck f 2 (n) (25) 

J \ / J —CO J Kmax 

20 



\f T R ff f = ldnp{n)( d ^ \ , (24) 



and now has four adjustable parameters. The last two terms constitute an endpoint fixing 
regulator. Without these terms, the smoothing regulator favors a linear PPS on the far left 
and the far right where the PPS is almost unconstrained by the CMB data. This behavior 
leads to spurious deviations from the fiducial model in the reconstructions. 

The reconstruction process may be characterized as follows. Let /j n be the input fractional 
variation of the PPS containing a number of features to be reconstructed. The recovered 
fractional variation is 

f recov f recov 

which is decomposed into a mean component f re cov and a random fluctuating component 
(of vanishing mean) 5f, whose covariance is given by 

(5f Sf) = ( J + R(X, a)) _1 7 (/ + R(X, a)) (27) 

and its mean is 

frecov = [Iff + Rff(\ a)]" 1 ////^ = Af in (28) 

The reconstruction operator A is a smoothing operator, which diminishes noise in f in from 
cosmic variance and detector fluctuations by averaging or smoothing, but also introduces a 
bias into the reconstruction f recov . This bias is unavoidable since some degree of smoothing 
is needed to avoid overfitting the data. 

For a numerical implementation, we must discretize f(k) by restricting to a finite- 
dimensional function space. We use a function space sufficiently large so that the results are 
insensitive to the dimension used. In the interest of numerical stability, it is advantageous 
to reduce the dimensionality of the function space as much as possible. We found that cubic 
5-splines are well suited for avoiding discretization artifacts while keeping the number of 
control points small. We expand as 

/(«) = X)/ i S i («). (29) 

i 

where B 1 {k) are cubic £>-spline basis functions. We work in the continuum limit where 
the roughness penalty suppresses variations on scales close to the control point spacing. 
Throughout a grid equally spaced in k = In A; with Ak = 0.01 is used. 
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(26) 



B. Reconstruction bias with other cosmological parameters fixed 



In this subsection we examine the mean reconstruction error when the other cosmological 
parameters (i.e. Hq, uj,, cj c , t, etc.) are fixed. In the next subsection we generalize taking into 
account uncertainties in these parameters and the following section examines the random 



error and questions of statistical significance. Fig. 10 plots the impulse (i.e., ^-function 
input) response of A without marginalizing over the cosmological parameters. The action 
of the reconstruction operator A on an impulse spreads a sharp localized input feature. 
Mathematically this is due to the smoothing regulator, which smears the variations in the 
power spectrum where its second derivative is large. 

The reconstruction operator introduces a bias in the reconstructed signal that under- 
estimates the height of features. Since the process is linear, a full description of the bias 
takes the form of the two-dimensional kernel, of which some representative cross sections are 



shown in Fig. 10 However two quantities provide a succinct partial characterization of the 
kernel: the impulse response width, defined as the standard deviation or square root of the 
second moment about the mean as a function of k in ; and the area under the curve, whose 
deviation from one indicates how much the amplitude of broad features is distorted. 



The limits ki ow and Khigh of the endpoint fixing regulator in eqn. (25), beyond which 
the spectrum is effectively constrained to coincide with the fiducial model, single out an 
interval of k space where to search for features. A concern is that deviations from the 
fiducial spectrum beyond these limits (where we cannot measure the power spectrum) could 
result in spurious features inside the free interval. We found that this not to be a problem 



when the other cosmological parameters were fixed. Several of the plots in Fig. [10] show the 
impulse response where the input is located in the fixed regions and we observe that except 
for some small artifacts near the boundary of a small integrated area, there are no artifacts 
in the interior of the free region. 



C. Reconstruction bias with uncertainty in cosmological parameters 

The situation studied in the previous subsection where the other cosmological parameters 
not associated with the PPS are assumed known with absolute precision is idealized. In 
practice, whether these parameters are determined internally from the CMB data alone or 

22 



k =0.002 Mpc" 1 



10" 2 10" 1 

k (Mpc' 1 ) 



k =0.4 Mpc" 1 




5.3 xlO" 
6.3 xlO" 
7.1 xlO" 



V./A ■ 



10" 2 10 

k (Mpc' 1 ) 



k =0.003 Mpc' 1 



/ 
/ 
/ 

/ 




\ 


0.04 

0.05 

0.07 


"■I ■ 


1 






V/ 


10" 2 10" 1 

k (Mpc' 1 ) 

fc c =0.01 Mpc' 1 



A 




0.97 

1.00 

1.01 








T^fvy. i i 




10" 2 10" 1 

k (Mpc' 1 ) 

fc c = 0.1 Mpc' 1 


i 


i 


1.05 

1.06 

1.00 


1 


L 


i 




i ^ 







i 


i 


0.11 

0.16 

0.21 




i 




1 1 1 1 1 L 





10"^ 10" 

(Mpc' 1 ) 

A; c =0.05 Mpc" 1 



10" 



k (Mpc' 1 ) 
fc c =0.3 Mpc" 1 



10" 1 




2.2 xlO 
2.7x10 
3.0 xlO 



10" z 10" 

k (Mpc' 1 ) 



1 


J 


i 





3.98/ 
3.99/ 
1.00/ 

// 










10" 2 10" 1 

(Mpc" 1 ) 




FIG. 10. Impulse response of reconstruction operator A (case with other cosmological 
parameters fixed). The curves correspond to A = 10 4 (blue), A = 10 3 (green) and A = 10 2 (red). 
In the shaded regions on the far left and right (where k < K[ ow = —5.5 or k, > Khigh = —1-5) an 
endpoint fixing prior (with a = 10 4 ) has been applied, pushing the reconstruction / to zero there. 
The numbers in the boxes indicate the area under the curves of the reconstruction kernel, which 
ideally should be as close as possible to one, the area under the delta function input. The values 
in the boxes in the two last plots show that there is little leakage from variations in the shaded 
regions. The lower right plot shows the impulse response width as a function of feature position. 
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whether additional non-CMB data is also used (acting as a prior on these parameters), some 
uncertainty in these additional parameters is always present. This uncertainty degrades the 
quality of the power spectrum reconstruction, because variations in the power spectrum can 
be traded against changes in these parameters. This subsection describes the properties of 
the mean reconstruction in the presence of these uncertainties, with and without external 
information. The next section treats the noise properties of the reconstruction. 

We found that for the case where the other cosmological parameters were determined from 
the CMB data alone, in order to prevent possible features in the fixed region from introducing 
spurious features in the free region and spurious shifts in the cosmological parameters, it is 
necessary to choose a wide interval [k,i ow , Khigh) = [10~ 4 , 0.3] Mpc -1 for the free region. ki ow 
here corresponds roughly to the current horizon scale (ki ow ~ l/f]o where 770 is the distance to 
the last scattering surface). k high is slightly above the k corresponding to £ ma x- The situation 
contrasts with the case of the previous subsection where the cosmological parameters are 
fixed, for which [ki ow , K high ] can be chosen freely without introducing artifacts. For the 
remainder of the paper we set a = 10 4 and define the free region as [10~ 4 , 0.3] Mpc -1 . 



The top row of Fig. 11 shows that the reconstruction operator badly distorts an infinitely 
narrow impulse input in the case where h, ouj,, and u c are uncertain and constrained only by 
the CMB data. This distortion is due to trading variations in the cosmological parameters 
against variations in the power spectrum, leading to a reconstruction kernel having long 
oscillating tails in the absence of additional non-CMB information, as in the case plotted 



here. Fig. 11 also shows the response to a Gaussian input profile of varying width as plotted 



in the bottom two rows. The bad behavior in the tails of the reconstruction kernel almost 



disappears when the kernel is applied to a broadened input. Fig. 12 shows that in the 
region of fc-space where the bias is most extreme, the response of an impulse away from the 
central peak oscillates rapidly as the location k c of the input impulse is varied. Although 
the amplitude of the response in the tail oscillates, its shape remains roughly constant as 
shown in the figure. This is why the oscillations in the tail interfere destructively almost 
canceling each other when the impulse input is replaced with a broadened Gaussian input. 

Allowing for uncertainty in h, u c and Ub, not only modifies the reconstruction kernel, but 
also causes features in the power spectrum to displace the maximum likelihood cosmological 
parameters. In the absence of a regulator the mean of this effect would vanish, but the 
mean is nonzero for A > 0. By expanding the transfer function to first order around the 
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FIG. 11. Smoothing Response to impulse and Gaussian inputs (with marginalization over 
the other cosmological parameters). The values of the prior parameters and their color mapping 

10~ 4 Mpc" 1 and 



are the same as in Fig. 
k high = 0.3 Mpc" 1 . 



10 The endpoint fixing regulator limits here are k, 



low 



input cosmological parameters, we can express the change in the CVs as a linear function of 
the change in the PPS and the change in the cosmological parameters 6(3: 



6C, = T cf f + T c ?6(3. 



(30) 



If we treat (/, 6(3) as a single parameter vector, we can write an extended reconstruction 
operator A that acts on the space of (/, 6(3), which we can write in block form as 



.4 



Aff A fP 
A/3f 



(31) 
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FIG. 12. Rapid oscillations in the marginalized smoothing operator A(k out , ki n ). The 
bottom plot shows A(k out , ki n ) versus k ou t for different ki n . The black vertical bars on the right 
indicate five values of ki n for which A(k out , ki n ) is plotted. We observe that the shape of the tail 
of the kernel remains approximately fixed but its overall amplitude undergoes oscillations as ki n 
is varied. The top figure shows A(k ou t, ki n ) as a function of ki n with k ou t fixed at the value of k 
corresponding to the position of the single vertical line on the left in the bottom plot. The points 
show the ki n values for the curves of the bottom plot. 



A is defined by inserting / from eqns. (22) or (23) into eqn. (28). Due to the off-diagonal 



block Apf, a feature in the PPS will in general shift the values of the cosmological parameters. 

For small A it is useful to include other data to prevent the other parameters from taking 
unreasonable values. The constraints from external data used here are h = 0.72 ± 0.08 from 
[36] (HST key project) and Q b h 2 = (2.2 ±0.1) x 10 -2 from [37J [38] (measured deuterium 
abundances), but Q c h 2 is left unconstrained by the external data. These constraints are 



implemented as Gaussian priors in our calculations. Fig. 13 shows that in the absence 
of external data, for small values of A the determination of the cosmological parameters 
degrades considerably, contradicting what is known from external data. With external data 
constraining h and w& the situation is remedied substantially, although for interesting values 
of A, up to 20% errors in Q c h 2 persist. 

With external datasets included, the deformation of a signal due to the reconstruction 



operator is lessened. In Fig. 14, the impulse and Gaussian response of the reconstruction 
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FIG. 13. Fractional error in cosmological parameters as a function of the smoothing 
parameter A. The solid curves show the errors with Gaussian priors on h and Qbh 2 , while 
the dashed curves show the errors without priors. At extremely large values of the smoothing 
parameter (A > 10 12 ), the form of the feature is restricted to a straight line, reducing the effective 
degrees of freedom in f(k) to two (a magnitude and tilt). In this case the error on the cosmological 
parameters is at its absolute minimum. At the other extreme, with A 55 10 _1 the form of the PPS 
is almost totally unconstrained and the constraints on the cosmological parameters are those from 
the external data. 



operator is shown. Compared to Fig. 11, the response is better behaved although still not 



as clean as the response when the cosmological parameters have been fixed (Fig. 10). This 



improvement is due to the decreased uncertainty in the cosmological parameters, which in 
turn lessens the need for the PPS to compensate for a shift in the cosmological parameters. 



Figs. [TT] and [14] demonstrates that features that are too narrow cannot be reconstructed 
faithfully because of confusion with variations of the other cosmological parameters. We de- 
fine the minimum reconstructible width AK mrw to describe this situation quantitatively. The 
quality of a reconstruction can be characterized by the mean square error, and setting the 
somewhat arbitrary threshold here to b = 10~ 2 , we deem a reconstruction to be satisfactory 
when 



J dK\Af Kc>a (K) - f Kc ,a(n)\ 2 < b J dK\f Kc>a (K) 



(32) 



Here f K a is a Gaussian centered at k c with standard deviation a. For a given k c , the 
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FIG. 14. Smoothing response to impulse and Gaussian inputs (with marginalization and 
external priors over the other cosmological parameters). The values of the prior parameters and 
their color mapping are the same as in Fig. 10 The endpoint fixing regulator limits here are 
how = 10" 4 Mpc" 1 and k high = 0.3 Mpc" 1 . 



minimum reconstructible width AK> mrw is denned as the smallest a such that the above 
condition still holds. 

Fig. [15] plots the minimum reconstructible widths for the marginalized and unmarginal- 
ized cases, including the impulse response width for the unmarginalized case for comparison. 
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FIG. 15. Minimum reconstructible width as defined in eqn. (32) with other parameters 
fixed (blue) and with marginalization including HST and deuterium external data (red). The 
impulse response width without marginalization has been plotted in green for comparison with the 
minimum reconstructible width. 
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IV. RECONSTRUCTION NOISE AND STATISTICAL SIGNIFICANCE 



The previous section discussed the mean (nonrandom) reconstruction error for features 
in the power spectrum superimposed onto a power law fiducial spectrum. In this section we 
examine the random error (having zero mean) of the reconstruction resulting from cosmic 
variance and other noise in the determination of the CMB angular power spectrum. We 
calculate the noise assuming the null hypothesis that the fiducial power spectrum is correct. 
Under this hypothesis, the noise covariance matrix is given by the expression 

N recon ={l + R(X, a)) (/ + R(X, a)) ~' (33) 

in the presence of the regulator or smoothing penalty. In the absence of a smoothing 
regulator, the noise would simply equal 7 -1 . If we had taken the Bayesian prior in the 
regulator seriously as a prior over the real inputs, the random error would be ( I+R(X, a) 



which is larger than eqn. (33 ). The difference between the two expressions can be understood 
as a consequence of the smaller than one eigenvalues of the operator + R(X, a)j I. 
The corresponding eigenvectors or eigenfunctions correspond to modes suppressed by the 
smoothness penalty. Since we are primarily interested in assessing the statistical significance 
of a potential first detection against the null hypothesis of a pure power law spectrum, we 



adopt the error in eqn. (33) in the calculations below. The operators above include both 
indices for the continuous parameter of / and the other cosmological parameters a in the 
case where these are not fixed. We report on both the random errors of the cosmological 
parameters and of the feature reconstruction. 
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A. Random error with other cosmological parameters fixed 
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FIG. 16. Random error with other cosmological parameters fixed. The impulse response 
width and random error autocorrelation function averaged over k G [0.01,0.1] Mpc" 1 is shown as 
a function of the smoothing parameter A. 



Fig. [16] shows how the impulse response width and reconstruction error (averaged over 
the k range [0.01,0.1] Mpc' 1 ) vary as a function of A. The bias decreases as A is lowered, 
but this entails a rising error from reconstruction noise, whose autocorrelation function is 



given by the diagonal elements N ream (k, k) as defined in eqn. (33). It can be seen from the 



plot that aAn irw is approximately constant. Fig. 17 shows the shape of the random error 
autocorrelation function. 

The reconstruction method was tested by generating several thousand Gaussian realiza- 
tions of a fiducial Ce power spectrum calculated using CAMB (35] for a given PPS and 
cosmological parameter values. The features had the Gaussian form 



f(k) = Aexp 



(In k — In k c 
2^ 



(34) 
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FIG. 17. Random error autocorrelation with other parameters fixed. The predicted 
error is dashed while the measured error is solid. 



parameterized by the width a, height A, and center k c , chosen according to A so that 
the feature is just wide enough to be accurately recovered (as determined by the impulse 



response width in Fig. 15) and just high enough to obtain a 5a deviation from the scale 



invariant spectrum. According to Figs. 15 and 17 the bias and error are minimal around 
O.OlMpc -1 < k < O.lMpc' 1 , so we test the reconstruction with features localized at k c = 
0.05 Mpc" 1 . 



Fig. 18 shows the result of 10 4 Monte Carlo simulations of the reconstruction of 4 mock 



features. The expected error in eqn. (33) agrees well with the standard deviation of the MC 



reconstructions in Fig. 17 
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FIG. 18. Reconstructed PPS from 10 4 Monte Carlo realizations with other parameters 
fixed for four different mock features. The actual PPS (solid red curve) and reconstructions 
average (dashed blue curve) agree. 
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B. Random errors with other cosmological parameters unfixed 
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FIG. 19. Minimum reconstructible width and random error autocorrelation averaged 
over k 6 [0.01, 0.1] Mpc" 1 with other parameters uncertain but with external priors on 
Hq and ujfj as function of A. 



In the case with uncertainty in the other cosmological parameters, the reconstruction 
becomes slightly more difficult. As in the fixed case, we reconstruct the PPS for four different 
features using 10 4 Gaussian random realizations of a fiducial Ci power spectrum, as shown 



in Fig. 20 The shape and position of the feature were chosen as in the fixed case except that 
the minimal reconstructible width is used instead of the impulse response width for choosing 
the input feature width. Comparing the average of all 10 4 reconstructed spectra (dashed 
blue line) and the actual PPS (solid red line), we see that the PPS is well reconstructed. 
Compared to the fixed case (Fig. |l8]) , the errors are slightly larger, particularly for k between 
O.OlMpc -1 and O.lMpc -1 . The reconstructed cosmological parameters with errors are given 



in Table 21 , both with and without external data Hq and 
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FIG. 20. Reconstruction of the PPS from 10 4 Gaussian random realizations of the Cg 
power spectrum for four different PPS. The actual PPS is shown as a solid red curve while 
the average of the reconstructions is shown as a dashed blue curve. The cosmological parameters 
have been marginalized over with the HST and deuterium external data included. 
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FIG. 21. Effect of external priors on Hq and cu on the cosmological parameters. Here 
for the external data we used the HST and deuterium data error bars but shifted the central values 
to coincide with the WMAP 7- year best fit values. It is the shift in these parameters as A is varied 
that is of interest here. 



Table 21 shows that without external data to constrain uj^ and H (from deuterium 
abundances and the HST key project results), the parameters take unreasonable and possibly 
negative values for small A. 
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C. Graphical representation of reconstruction 



In this subsection we present a graphical representation of the reconstructed power spec- 



trum with errors as given in eqn. (33). As we saw, f(k) is not measurable at a point. 



Rather it is the integral over the kernel that is being reconstructed. We therefore indicate 
error boxes whose width is equal to the impulse response width for the fixed case and the 
minimal reconstructible width for the case of uncertainty in the other parameters. 



Figs. 22 and 23 illustrates the reconstruction of the PPS without the other parameters 
fixed and with the other parameters uncertain but with the external data included, respec- 



tively. The input feature shapes are the same as in Figs. 18 and 20 The la and 2a errors 
in the reconstruction are represented by the height of the green shaded boxes. Thus the 
PPS has been successfully reconstructed if the actual PPS passes through most of the error 



boxes, which is the case in Figs. 22 and 23 
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FIG. 22. Reconstruction of mock features with other cosmological parameters fixed. 

The actual feature (solid red) agrees with the average reconstruction (dashed blue). Green boxes 
indicate the error of the reconstruction. Their heights represent the la and 2a errors and their 
widths represent the impulse response width (as evaluated at the box center). 
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FIG. 23. Reconstruction with h, LJb, and uj c uncertain with HST and deuterium abun- 
dance external constraints included. 



37 



D. Non-Gaussian corrections 

The preceding discussion assumed a quadratic likelihood. However the actual likelihood 
is nonlinear mainly because of non-linearity in the \ow-£ likelihood for the C/s and the non- 
linear dependence of the predicted C/s as a function of the cosmological parameters. Once 
the cosmological parameters have been fixed, the derivative of the likelihood with respect 
to variations in the power spectrum is easy to compute analytically. Moreover, many of 
the high-£ likelihoods for a theoretical Ci spectrum taking into account partial sky cover- 
age, nonuniform instrument noise, and removal of high-£ residual galactic and extragalactic 
contamination are based on approximations for which derivatives of the likelihood with re- 
spect to changes in the C/s may be calculated with simple code modifications. For feature 
reconstruction it is the high-£ likelihood, and not the \ow-£ likelihood, where a pixel-based 
approach is required, that is most relevant. Consequently, for obtaining the maximum likeli- 
hood solution, it is feasible and simple to apply the Newton-Raphson root finding algorithm 
to obtain the exact ML reconstruction within a few iterations despite the massive increase 
in the number of parameters. Accounting for non-Gaussianity in the fluctuations about the 
ML reconstruction is more difficult. MCMC methods are not feasible because there are too 
many parameters. There is likely to be mild non-Gaussianity impacting on the analysis in 
Section 4. This is an issue currently under investigation. 
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V. CONCLUSION 



In this paper we have explored a model-independent scheme for reconstructing possible 
features in the primordial power spectrum based on a penalized likelihood. After exploring 
the structure of the available information, we study the properties of the reconstruction 
for a range of degrees of smoothing. What features in the primordial power spectrum 
can be reconstructed and with what statistical significance is quantified for a collection of 
representative shapes. We found that very narrow features cannot be reliably reconstructed 
due to confusion with the cosmological parameters even when external data is used to fix H 
and Ub- It will be interesting to see whether future data will uncover significant deviations 
from a power law adiabatic primordial power spectrum. 
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